#include <iostream>
#include <cstdlib>
#include <cmath>
#include <iomanip>
#include <limits>
#include "parse.hpp"
#include "rat.hpp"

typedef double(*double_fun)(const double&); // define type "double_fun"
typedef long double myfloat;
double trapz1(const double_fun& f, 
	      const double& a, const double& b, 
	      const size_t N){
	double h = (b-a)/N;
	double sum = (f(a)+f(b))/2;
	for(size_t i=1; i<N; i++){
		sum += f(a+i*h);
	}
	return h*sum;
}
template <typename Function>
double trapz2(const Function& f, 
	      const double& a, const double& b, 
	      const size_t N){
	double h = (b-a)/N;
	double sum = (f(a)+f(b))/2;
	for(size_t i=1; i<N; i++){
		sum += f(a+i*h);
	}
	return h*sum;
}

template <typename T, typename Function>
T trapz3(const Function& f, 
	      const T& a, const T& b, 
	      const size_t N){
	T h = (b-a)/N;
	T sum = (f(a)+f(b))/2;
	for(size_t i=1; i<N; i++){
		sum += f(a+i*h);
	}
	return h*sum;
}

template <typename T=double> // T defaults to type double
struct power
{
	T m;
	power(const T m=2) : m(m) { }
	T operator()(const T& x) const
	{
		using std::pow;
		return pow(x, m);
	}
};
double x2(const double& x)
{
	return x*x;
}

auto square_fun = [](const double& x) { return x*x; };

int main(int argc, char *argv[]) {
	long double exact = 1.0L/3.0L;
	size_t N = 128;
	if(argc > 1) parse(argv[1], N);
	std::cout
		<< std::setprecision(std::numeric_limits<long double>::digits10+1)
		<< std::scientific;
	myfloat Q = trapz3<myfloat>(power<myfloat>(), 0, 1, N);
	using std::abs;
	std::cout
		<< N << "\t"
		<< Q << "\t"
		<< abs(Q - exact) // absolute error
		<< std::endl;

}



